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1.  Technical  Objective 

This  research  is  considering  the  technologies  needed  to  realize  the  promise  of  optimum 
rates  of  convergence  of  adaptive  high  order  methods  based  on  hp-discretizations.  The 
focus  of  the  current  project  was  the  development  of  the  technologies  needed  to  support 
robust  p-version  mesh  generation  of  general  non-manifold  geometric  domains  of  interest 
to  the  Navy. 

2.  Technical  Approach 

To  meet  the  goal  of  providing  appropriate  meshes  for  p-version  finite  element  analysis  for 
general  domains,  an  automatic  mesh  generation  procedure  to  create  curved  finite  element 
meshes  with  sufficient  order  of  geometric  approximation  based  on  the  order  of  polynomial 
used  in  the  finite  element  basis  was  created.  Although  the  initial  straight-sided  mesh  gen¬ 
erated  by  automatic  mesh  generators  has  all  valid  elements,  the  process  of  curving  the 
mesh  edges  and  faces  to  the  appropriate  model  boundaries  often  yields  elements  with 
invalid  shapes.  Therefore,  the  developed  procedures  perform  mesh  modifications  to  pro¬ 
duce  a  set  of  valid  curved  elements. 

A  key  ingredient  is  the  geometric  representation  of  the  mesh  entities.  The  standard  method 
of  Lagrangian  interpolation  does  not  lend  itself  to  effective  procedures  as  the  order  of  the 
elements  increases.  Therefore,  an  alternative  geometric  form  based  on  Bezier  approximat¬ 
ing  geometry  was  developed. 

Many  of  the  physical  domains  of  interest  to  the  Navy  are  structures  dominated  by  thin  sec¬ 
tions  for  which  the  desired  finite  element  discretizations  employ  thin  volume  elements. 
For  the  most  effective  application  of  p-version  finite  elements  it  is  desirable  to  generate 
elements  in  these  thin  volumes  that  do  not  have  long  diagonals  going  through  the  “thick¬ 
ness  direction”.  This  creates  problems  for  current  automatic  mesh  generators  that  use  only 
tetrahedra.  Therefore,  procedures  were  developed  to  eliminate,  to  the  greatest  possible 
extent,  such  diagonals. 

3.  Results 

3.1  Procedure  for  the  generation  of  curved  meshes 

The  procedure  to  curve  the  elements  of  an  initially  straight-sided  mesh  employ  a  set  of 
mesh  curving  tools  based  on  standard  operators  and  a  set  of  mesh  modification  operators 
that  are  capable  of  modifying  the  local  mesh  topology.  The  full  set  of  mesh  modification 
operators  includes  swap,  split,  and  collapse  operators  that  can  be  applied  to  various  mesh 
entities.  In  addition  to  those  applied  as  independent  operations,  there  are  times  when  they 
can  be  applied  in  various  combinations.  Finally,  there  is  a  general  cavity  creation  and  tri¬ 
angulation  procedure  that  can  be  applied  as  needed. 
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These  procedures  [6,10]  have  been  implemented  in  the  mesh  generation  tool  kit  [7]  and 
have  been  shown  to  provide  superior  results  to  previous  more  ad-hoc  approaches  [2,3]. 
Experience  on  the  geometry  correction  of  straight-sided  meshes  (that  is,  the  snapping  of 
refinement  points  [4])  indicates  that  careful  analysis  of  the  local  situation  is  most  effective 
in  determining  which  mesh  modification  operations  to  apply  to  ensure  there  is  progress  in 
improving  the  mesh.  The  development  of  these  procedure  to  curved  elements  is  compli¬ 
cated  by  the  need  to  do  the  checks  against  the  curved  element  geometry  which  does  intro¬ 
duce  substantial  computational  cost.  Figure  1.  shows  the  application  of  this  procedure 
where  both  quadratic  and  cubic  element  geometries  are  used,  while  Figure  2.  shows  the 
application  of  the  procedure  on  a  mechanical  part. 


Model  Geometry 


Original  Linear  Mesh 


Figure  1.  Mesh  curving  example  using  quadratic  and  cubic  element  geometry 


Figure  2.  Mesh  curving  for  a  mechanical  part 


3.2  Geometric  Representation  of  Curved  Mesh  Entities 

The  two  basic  approaches  to  represent  the  geometry  of  mesh  edges  and  faces  on  the  model 
boundary  are: 

•  Exact  model  geometry  through  geometric  model  interrogation.  This  method  has  been 
previously  implemented  and  found  to  be  computationally  expensive  [1]. 


•  A  geometric  approximation  of  sufficient  order  to  ensure  proper  convergence  of  the 
finite  element  solution. 

For  mesh  edges  on  the  interior  of  the  model  that  need  to  be  curved  to  preserve  the  shape  of 
attached  already  curved  elements,  any  polynomial  fit  of  an  order  consistent  with  the  finite 
element  basis  of  the  element  is  satisfactory. 

Two  options  for  the  geometric  approximation  of  curved  boundary  and  interior  mesh  edges 
and  faces  considered  are  Lagrange  interpolation  and  Bezier  approximating  polynomials. 
Although  the  definition  of  Lagrange  polynomials  is  straight  forward,  the  determination  of 
interference  and  control  of  these  polynomials  for  higher  than  quadratic  geometry  becomes 
difficult.  Therefore,  efforts  on  this  project  focused  on  the  development  of  the  Bezier  geo¬ 
metric  approximation  [5,6,9].  The  key  Bezier  geometry  properties  of  interest  are: 

•  Can  be  as  high  a  degree  as  desired 

•  Convex  hull  provides  smoother  and  more  controllable  approximation 

•  Better  properties  to  allow  more  efficient  intersection  checks 

•  Derivatives  and  products  of  Beziers  are  also  Beziers 

•  Efficient  algorithms  for  degree  elevation  and  subdivision 

The  key  technical  issues  associated  with  the  implementation  of  Bezier  element  geometry 
are: 

•  Defining  mesh  entity  shapes  on  the  model  boundary  based  on  interactions  with  the  geo¬ 
metric  modeling  system.  This  includes  determining  the  points  on  the  model  entity  that 
the  mesh  entity  should  pass  through,  the  corresponding  parametric  locations  on  the 
mesh  entity  for  each  interpolation  point,  and  if  the  resulting  shape  self-intersects  and 
how  to  correct  the  shape  if  needed 

•  Setting  mesh  entity  order  based  on  entity  closure  order  requirement 

•  Ensuring  the  validity  of  the  elements 

•  Geometric  properties  and  operations  to  use  in  the  process  of  correcting  mesh  entity 
shapes 

•  Providing  an  interface  for  the  mesh  curving  routine 

The  mesh  entity  Bezier  approximations  must  be  constructed  from  pointwise  interrogations 
of  the  geometric  model.  The  process  must  account  for  issues  associated  with  geometric 
modeling  systems  face  parametric  coordinates  including  periodic  faces  and  edges,  degen¬ 
erate  points  on  faces  and  parametric  system  distortion.  Mesh  edges  classified  on  model 
face  are  parameterized  based  on  a  cord  length  procedure  to  control  the  parametric  space 
distortion.  Procedures  to  deal  with  parametric  space  periodicity  and  typical  degeneracies 
have  been  defined. 

The  next  key  task  was  determination  and  elimination  of  any  self-intersections  that  arise. 
Self-intersections  can  occur  in  the  interpolation  geometry  when  the  model  geometry  is  of 
higher  degree  than  the  degree  of  the  mesh  entity  or  is  a  piece-wise  shape  (such  as  a 
NURB).  Inexpensive  methods  for  determining  potential  self-intersections  are  being 
devised.  In  the  case  of  mesh  edges,  the  method  makes  use  of  the  variation  diminishing 
property  of  Bezier  curves,  which  states  that  an  infinite  plane  can  not  intersect  a  Bezier 
curve  more  times  than  it  intersects  its  control  polyline.  This  allows  the  development  of  a 
procedure  to  find  obvious  edge  intersection  problems. 

Detecting  self-intersections  in  mesh  faces  is  more  difficult  since  Bezier  Surfaces  do  not 
exhibit  the  variation  diminishing  property.  The  process  begins  with  an  analogous  test 
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defined  for  Bezier  Surfaces  based  on  the  projection  of  the  control  polygon  onto  the  planar 
face  of  the  linear  mesh  face.  If  the  projection  of  the  control  polygon’s  points  form  a  non¬ 
intersecting  tessellation  the  surface  is  less  likely  to  self-intersect.  This  builds  on  the  fact 
that  the  control  polygon  resulting  from  continually  degree  elevating  a  surface  will  in  the 
limit  become  the  Bezier  Surface. 

The  previous  implementations  forced  all  mesh  entities  to  be  raised  to  the  highest  degree 
requested.  Since  many  models  have  planar  faces  and  linear  edges,  this  restriction  results  in 
mesh  entity  shapes  that  are  more  complex  than  needed.  By  allowing  mesh  entities  to  have 
different  polynomial  degrees,  mesh  entity  shapes  are  as  complex  as  required.  In  order  to 
allow  mesh  entities  to  have  different  polynomial  degrees,  a  mesh  entity  may  need  to  pro¬ 
vide  higher  order  representations  in  order  for  another  mesh  entity  that  has  the  entity  as  part 
of  its  boundary  to  be  of  higher  order. 

Determining  the  validity  of  region  elements  requires  that  the  determinate  of  the  Jacobian 
be  positive.  Properties  of  Bezier  regions  useful  in  this  process  are: 

•  Partial  derivatives  of  the  region  are  themselves  Bezier  functions 

•  The  Jacobian  determinant  is  defined  by  box-product  of  partial  derivatives 

•  Since  the  product  of  2  Bezier  functions  is  also  a  Bezier  function,  the  Jacobian  determi¬ 
nant  is  also  a  Bezier  function 

•  In  the  case  of  a  tetrahedron  the  function  is  of  order  3(p-l),  where  p  is  the  order  of  the 
original  shape 

•  Since  a  Bezier  function  is  bounded  by  its  convex  hull,  the  Jacobian  determinant  func¬ 
tion  inside  the  region  is  bounded  by  the  convex  hull  of  its  control  points  (which  in  this 
case  are  scalars) 

•  A  region  is  valid  globally  if  the  minimum  control  point  of  the  Jacobian  determinant 
function  is  positive 

The  box  product  terms  that  compose  the  Jacobian  determinant  function  can  be  used  to 
determine  how  a  region  should  be  corrected.  An  outline  of  the  procedure  is 

•  For  each  Jacobian  that  is  negative: 

•  Identify  the  minimum  Box  Product  term  contributing  to  the  Jacobian 

•  Identify  the  region  control  points  involved  in  the  box  product  vectors  that  can  be 
moved.  Control  Points  may  be  constrained  due  to  being  associated  with  mesh  entities 
on  the  model  boundary  or  constrained  to  prevent  other  mesh  regions  from  becoming 
invalid 

•  Identify  the  minimum  angle  to  make  the  Box  Product  Term  positive 

•  Determine  which  control  points  of  region  that  defines  the  vector  should  be  displaced  in 
order  to  rotate  the  vector 

•  If  this  change  would  result  in  invalidating  a  neighboring  region  it  must  be  modified  to 
accommodate  the  shape  change 

•  If  the  region  is  still  invalid  then  perform  one  of  the  following: 

-  Degree  elevate  the  region’s  shape  if  possible  in  order  to  increase  the  degrees  of 
freedom  available 

-  Sub-divide  the  region  in  order  to  refine  the  mesh  and  introduce  more  degrees  of 
freedom 

After  the  surface  mesh  has  been  properly  curved  and  all  the  mesh  regions  have  been  made 
valid,  the  last  step  in  the  mesh  curving  process  is  to  detect  global  mesh  intersections.  This 
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procedure  uses  the  convex  hull  property  of  a  Bezier  curve,  surface,  or  volume  which  states 
it  is  contained  in  the  convex  hull  formed  by  its  control  points.  Therefore,  if  the  convex 
hulls  of  two  mesh  entities  do  not  intersect  then  the  entities  do  not  intersect.  If  the  convex 
hulls  do  intersect  then  the  convex  hulls  can  be  refined  by  subdivision.  Refinement  contin¬ 
ues  until  either  it  is  determined  that  the  hulls  do  not  intersect  or  that  an  intersection  has 
occurred  (within  an  acceptable  tolerance).  Once  an  intersection  is  detected,  the  control 
points  of  the  interfering  geometry  are  displaced  in  order  to  remove  the  intersection. 

To  isolate  the  details  of  the  Bezier  mesh  geometry  representation  from  the  mesh  curving 
and  modification  procedures,  all  interactions  with  the  shape  information  is  through  general 
classes.  The  use  of  the  topologically  specific  sub-classes  allows  for  the  effective  account¬ 
ing  for  the  local  parametric  spaces  of  the  mesh  entities.  Every  subclass  shares  the  same 
interface  with  degree,  mesh  entity  and  a  set  of  control  points  location.  The  information  in 
the  element  entity  shape  subclasses  is  operated  on  by  operators  which  are  used  in  the  pro¬ 
gram  logic  that  curves  the  mesh  entities  and  performs  the  appropriate  mesh  modification 
operations  to  ensure  the  final  curved  mesh  is  valid. 

In  addition  to  the  ability  to  create  curved  meshes,  a  number  of  issues  related  to  the  details 
of  the  geometric  approximation  used,  ranging  from  the  basic  order  used  to  the  selection  of 
interpolation  points  for  fitting  the  Bezier  geometry,  must  be  considered.  Although  theoret¬ 
ical  studies  clearly  show  a  loss  of  rate  of  convergence  when  the  geometric  approximation 
is  too  low,  many  argue  that  for  purposes  of  “engineering  accuracy”  only  quadratic  element 
geometry  is  needed.  The  simple  studies  presented  in  references  [5,6]  demonstrate  this  is 
not  true  and  that  very  poor  results  are  obtained  when  the  order  of  the  basis  exceeds  that  of 
the  element  geometric  approximation  for  coarse  meshes.  These  studies  also  indicate  the 
importance  of  accounting  for  the  curvature  for  selecting  the  distribution  of  interpolation 
points  for  fitting  the  Bezier  geometry. 

3.3  Specialized  Mesh  Generation  for  Thin  Sections 

Many  Naval  structures  are  constructed  as  an  assembly  of  stiffened  shells  where  the  thick¬ 
ness  of  the  individual  components  can  vary  through  the  structure.  One  approach  to  the  p- 
version  acoustics  analysis  of  these  structures  is  to  represent  the  thickness  as  a  volume 
component  in  the  model  and  use  volume  elements  where  the  p-value  in  the  thickness 
direction  is  low.  This  allows  the  direct  generation  of  a  mesh  of  volume  elements  in  both 
the  portion  of  the  water  domain  represented  and  in  the  thin  volumes  of  the  shell  structures. 
The  problem  with  this  approach  is  the  meshing  procedures  use  tetrahedra  elements  so 
there  would  be  long  through  the  thickness  diagonals  created  that  degrade  the  ability  to 
apply  effectively  p-version  representations  through  the  thickness.  An  alternative  approach 
is  to  represent  the  thin  sections  as  faces  in  the  geometric  model  with  thickness  attributes 
associated  with  them.  The  needed  meshes  in  the  thin  sections  would  then  be  generated  by 
first  getting  a  curved  surface  mesh  for  all  those  faces  (see  an  example  in  Figure  3.).  The 
volume  meshes  would  then  be  created  through  the  application  of  specialized  meshing 
operations  that  avoid  long  through  the  thickness  diagonals.  The  key  to  the  development  of 
this  approach  is  the  effective  treatment  of  the  face  intersections  where  more  than  two  faces 
share  the  edge. 

A  study  of  typical  stiffened  shell  intersection  was  undertaken  to  determine  the  most  appro¬ 
priate  mesh  configuration  which: 

•  Do  not  alter  the  desired  surface  mesh 
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Model  Mesh 

Figure  3.  Curved  non-manifold  mesh  of  domain  with  thin  sections. 


•  Create  the  minimal  number  of  volume  elements 

•  Do  not  have  long  diagonal  mesh  edges  in  the  through  thickness  direction 

The  alternative  approaches  to  perform  such  meshing  procedures  indicated  the  need  for  the 
following  mesh  topologies: 

•  hexahedron 

•  wedge 

•  pyramid 

•  five-noded  created  by  collapsing  one  of  the  triangular  faces  of  a  wedge  to  a  edge 

•  seven-noded  created  by  collapsing  one  edge  causing  two  of  the  faces  to  become  of  the 
triangular 

The  construction  of  the  element  shape  functions  for  all  these  element  topologies  is  to 
employ  the  shape  function  decomposition  approach  of  reference  [8].  Note  that  the  pyra¬ 
mid,  five-noded  degenerate  wedge  and  7-noded-degenerate  hexahedra  do  not  have  an 
appropriate  natural  coordinate  system  available.  In  each  case  the  element  must  be  con¬ 
structed  by  the  appropriate  degeneration  of  the  parametric  system  for  the  hexahedron, 
wedge  and  hexahedron  respectively.  This  introduces  a  second  mapping  that  must  be 
accounted  for  in  the  construction  of  the  region  blend  functions  and  in  performing  numeri¬ 
cal  integration.  The  region  blends  are  constructed  by  mapping  from  the  original  parent  ele¬ 
ment  to  the  degenerate  form,  defining  the  inverse  mapping  and  substituting  this  as 
appropriate  into  the  original  region  functions  for  the  non-degenerate  form.  Note  the  inte¬ 
gration  process  must  integrate  over  the  degenerate  form.  Therefore,  the  Jacobian  of  the 
mapping  from  the  non-degenerate  parametric  form  to  the  degenerate  form  appears  in  the 
integrand  evaluation. 

4.  References 

[1]  Dey,  S.,  Shephard,  M.S.  and  Flaherty,  J.E.,  “Geometry-based  issues  associated  with  p- 
version  finite  element  computations”,  Comp.  Meth.  Appl.  Meek  Engng.,  150:39-55, 
1997. 

[2]  Dey,  S.,  O’Bara,  R.M.  and  Shephard,  M.S.,  “Curvilinear  mesh  generation  in  3D”, 
Proc.  8th  Int.  Meshing  Roundtable,  SAND  99-2288,  pp.  407-417,  1999. 

[3]  Dey,  S.,  O’Bara,  R.M.  and  Shephard,  M.S.,  “Curvilinear  mesh  generation  in  3D”, 
Computer-Aided  Design,  33:199-209, 2001. 

[4]  Li,  X.,  Shephard,  M.S.  and  Beall,  M.W.,  “Accounting  for  curved  domains  in  mesh 
adaptation”,  to  appear  International  Journal  for  Numerical  Methods  in  Engineering, 
2003. 


6 


[5]  Luo,  X.,  Shephard,  M.S.  and  Remade,  J.-F.,  “The  Influence  of  Geometric  Approxima¬ 
tion  on  the  Accuracy  of  High  Order  Methods”,  Numerical  Grid  Generation  in  Compu¬ 
tational  Field  Simulations  2002,  Engineering  Research  Center  for  Computational 
Field  Simulations,  Mississippi  State  University,  Mississippi,  pp.  189-198,  2002. 

[6]  Luo,  X.,  Shephard,  M.S.,  Remade,  J.-F.,  O’Bara,  R.M.,  Beall,  M.W.,  Szabo,  B.A.  and 
Actis,  R.,  “p- Version  Mesh  Generation  Issues”,  11th  International  Meshing  Round¬ 
table,  Sandia  National  Laboratories,  pp.  343-354,  2002. 

[7]  Shephard,  M.S.,  “Meshing  environment  for  geometry-based  analysis”.  International 
Journal  for  Numerical  Methods  in  Engineering,  47(1-3):169-190,  2000. 

[8]  Shephard,  M.S.,  Dey,  S.  and  Flaherty,  J.E.,  “A  straight  forward  structure  to  construct 
shape  functions  for  variable  p-order  meshes”,  Comp.  Meth.  Appl  Mech.  Engng., 
147:209-233,  1997. 

[9]  Shephard,  M.S.,  O’Bara,  R.M.,  Luo,  X.  and  Li,  X.,  “Algorithm  to  generate  p-version 
meshes  for  curved  domains”,  p  and  hp  Finite  Element  Methods:  Mathematics  and 
Engineering  Practice,  Washington  University,  St.  Louis,  MO,  pp.  103-104,  2000. 

[10] Shephard,  M.S.,  Luo,  X.  and  Remade,  J.F.,  “Meshing  for  p-version  finite  element 
methods”,  WCCM  V:  Fifth  World  Congress  on  Computational  Mechanics,  Vienna 
Univ.  of  Tech,  Vienna,  Austria,  p.  1-464,  2002. 


7 


